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ABSTRACT 

We have recently introduced a novel statistical measure of dark matter clustering 
in phase space, the particle phase space average density ( P 2 SAD ). In a two-paper 
series, we studied the structure of P 2 SAD in the Milky-Way-size Aquarius haloes, 
constructed a physically motivated model to describe it, and illustrated its potential 
as a powerful tool to predict signals sensitive to the nanostructure of dark matter 
haloes. In this work, we report a remarkable universality of the clustering of dark 
matter in phase space as measured by P 2 SAD within the subhaloes of host haloes 
across different environments covering a range from dwarf-size to cluster-size haloes 
(10 10 — 10 15 Mq). Simulations show that the universality of P 2 SAD holds for more 
than 7 orders of magnitude, over a 2D phase space, covering over 3 orders of magnitude 
in distance/velocity, with a simple functional form that can be described by our model. 
Invoking the universality of P 2 SAD, we can accurately predict the non-linear power 
spectrum of dark matter at small scales all the way down to the decoupling mass limit 
of cold dark matter particles. As an application, we compute the subhalo boost to the 
annihilation of dark matter in a wide range of host halo masses. 
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1 INTRODUCTION 

The Cold Dark Matter (CDM) paradigm of structure forma¬ 
tion predicts the clustering of dark matter (DM) into gravi¬ 
tationally bound haloes in a very large range of scales, from 
the decoupling of CDM particles (e.g., ICC 11 — 10 -3 Mq, 
Bringmann 2009, depending on the model) to the massive 
10 15 Mq clusters virialized recently. Because in CDM most 
of the DM today is predicted to be inside haloes, accurately 
following the evolution of the DM phase space distribution 
within these highly non-linear structures is a challenging 
task. The tremendous progress of numerical A—body simu¬ 
lations have made it possible to cover the dynamical range 
paramount to galaxy formation, from large (~Mpc) to sub- 
galactic (~100 pc) scales, but is not yet feasible to explore 
the DM clustering at even lower scales, which we refer to as 
the nanostructure of DM haloes. This unresolved regime is 
however of prime importance in experiments searching for 
non-gravitational DM signatures that use theoretical pre¬ 
dictions, which in many cases are quite sensitive to the 
nanostructure of DM haloes. For instance, the DM self- 
annihilation rate over an entire halo is dominated by events 
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occurring in the plethora of its unresolved subhaloes, while 
the scattering rate in direct detection experiments is sensi¬ 
tive to the local DM phase space distribution. 

Even if there are no new DM interactions with visible 
particles, it might still be possible to detect the DM clus¬ 
tering at small scales through the gravitational influence of 
DM on astrophysical sources. Although haloes with a mass 
below the limit for atomic cooling are expected to be de¬ 
void of cold gas and stars, their presence can still be made 
evident through the gravitational lensing they produce on 
background sources. For instance, subhaloes close to the 
atomic cooling limit (~ 10® — 10 8 Mq) might be responsi¬ 
ble for the lensing flux-ratio anomalies observed in quasars 
(e.g. Xu et al. 2015), while even smaller subhaloes could be 
probed by time-varying lensing effects as they move along 
the line of sight of cosmological sources such as pulsars and 
quasars (Baghram et al. 2011; Rahvar et al. 2014). 

The spatial DM clustering at subhalo scales is known to 
have nearly universal properties within the range of scales 
resolved in simulations: a smooth radial distribution de¬ 
scribed by a spherically averaged two-parameter NFW pro¬ 
file (Navarro et al. 1997), and a hierarchy of subclumps (well 
described by NFW profile truncated at the tidal radius) with 
an abundance that is a power-law in mass, the subhalo mass 
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Simulation 

A /200 [Mq] 

^200 [Mpc] 

m p [Mq] 

e[kpc] 

A/s u b [Mg] 

A aub [xl0 5 ] 

Reference 

Dwarf-size 

1.3 x 10 10 

0.047 

1.2 x 10 3 

0.034 

3.2 x 10 s 

2.69 

Vogelsberger et al. (2014) 

Milky-Way-size 

1.8 x 10 12 

0.246 

1.4 x 10 4 

0.066 

1.3 x 10 11 

96.14 

Springel et al. (2008) 

Group-size 

8.7 x 10 13 

0.939 

1.5 x 10 7 

1.6 

6.8 x 10 12 

4.55 

this work 

Cluster-size 

2.9 x 10 15 

3.022 

2.4 x 10 8 

4.0 

1.5 x 10 14 

6.29 

this work 


Table 1. Summary of some of the properties of the simulations we analyze (at 2 = 0): the virial mass (M 200 ) and radius (r20o), the 
particle mass (m p ), the Plummer-equivalent maximum physical softening length (e), the total mass (M sll n)- the total number of particles 
(JV au b) in subhaloes resolved with more than 20 particles, and the reference to each simulation. Except for the Group-size case, all 
simulations have a lower resolution version (by a factor of 2 in softening), that we use to test convergence in P 2 SAD. 


function (e.g. Springel et al. 2008). Although the properties 
of the hierarchy of subhaloes are complicated by the time- 
varying effects of tidal disruption, (extrapolating) these ob¬ 
served universalities are the basis of most predictions of the 
DM clustering at small, unresolved scales. 

A complimentary picture emerges by considering the 
DM clustering in phase space, which becomes particularly 
relevant in potential signals of new DM interactions that are 
sensitive to the DM velocity distribution. In a series of pa¬ 
pers (Zavala & Afshordi 2014a,b, henceforth Papers I and 
II), we introduced the two-dimensional particle phase space 
average density ( P 2 SAD ), a coarse-grained phase space den¬ 
sity, which is a novel measure of DM clustering at small 
scales. In Paper I, we found signs of a (near-)universality 
in the structure of P 2 SAD over assembly history and red- 
shift for Milky-Way-size haloes, while in Paper II we pre¬ 
sented a physical model of P 2 SAD based on the stable 
clustering hypothesis in phase space (Davis & Peebles 1977; 
Afshordi et al. 2010), the spherical collapse model, and tidal 
disruption of subhaloes. We then showed how this model can 
be used to predict DM annihilation signals. 

In this work, we investigate further the structure of 
P 2 SAD and find that is remarkably universal for DM inside 
subhaloes across a wide range of host halo masses. It has a 
functional form that can be described by a simple paramet¬ 
ric formula with a structure that is motivated (modelled) by 
our physical prescription. Such universality in phase space 
makes it a powerful tool to describe the DM clustering at 
small unresolved scales, while its simplicity makes it useful 
to predict several potentially observable DM signals. 


2 PARTICLE PHASE SPACE AVERAGE 
DENSITY ( P 2 SAD ) IN SUBHALOES 

We follow the same notation as in Paper I to define 
P 2 SAD = 3(Ax, Au)v 6 as the mass-weighted average (over 
a volume V6 in phase space) of the coarse-grained phase 
space DM density, on spheres of radius Ax and An, in posi¬ 
tion and velocity spaces, respectively: 


3{Ax,Av) Vb 


J Vg d 3 xd 3 v/(x, v)/(x + Ax, v + Av) 
/ Ve d 3 xd 3 v/(x,v) 
</(x,v)/(x + Ax,v + Av )) V6 

</>v„ ( j 


where /(x, v) is the phase space distribution function at the 
phase space coordinates x and v, and (f) Ve is the average 


phase space density in the volume V6' 

fv e d 3 xd 3 v/(x, v) _ Mv(j 

">v e - y 6 - y 6 


(2) 


In a simulation, the dark matter density field is repre¬ 
sented by a discrete set of N particles, each with a mass m v . 
In this representation, Eq. (1) is estimated as: 

m p (N p ( Ax,Av)) Ve 
Ve{Ax,Av) 


^(Ax,Au) s im — 


(3) 


where (N p ) is the average number of simulation particles 
(over Ve) within shells of thickness ( 8x,8v ) at a radius Ax 
and Av in phase space centred on each of the particles, and 
Ve(Ax, Av) is the phase space volume of a given shell (see 
Fig. 1 of Paper I for a diagram illustrating how P 2 SAD is 
estimated in a simulation). 

In Papers I and II, we studied P 2 SAD averaged over 
all particles within the virial radius of each of the Aquarius 
haloes (Springel et al. 2008). We define the virial radius as 
the radius enclosing a sphere with mean density 200 times 
the critical value (r 2 oo). Here, we focus instead only on the 
particles inside self-bound subhaloes within the virial ra¬ 
dius of the host, and extend the analysis to different host 
halo masses, covering a range from dwarf-size to cluster-size 
haloes (10 10 — 10 15 Mg, see Table 1 for a summary of the 
simulations we use). 

We find that P 2 SAD is remarkably universal across all 
halo masses (see Fig. 1). In Paper I, we had already re¬ 
ported some level of universality for the limited mass range 
of the Aquarius haloes (~ 8 x 10 11 — 2 x 10 12 Mq), and 
for P 2 SAD averaged over all particles within r 2 oo- Here we 
show that the universal character of P 2 SAD within sub¬ 
haloes extends across the more than 5 orders of magnitude 
in halo masses that we have explored. Notice the large range 
of scales (3 orders of magnitude in both, distance and ve¬ 
locity separations; 18 orders of magnitude in the 6D phase 
space volume) where the universality of P 2 SAD holds, while 
varying by 7 orders in magnitude. 

Subhaloes with a characteristic physical scale and veloc¬ 
ity contribute more prominently to a fixed value of P 2 SAD, 
i.e., each contour in Fig. 1 is mainly representative of sub¬ 
haloes of similar size. The maximum scale that subhaloes 
can have is limited by the size of their host (governed by 
the tidal stripping in the parent halo), and thus, above 
this maximum scale, subhaloes of a given host cannot con¬ 
tribute to low values of P 2 SAD, which is why the subhaloes 
in the dwarf (Milky-Way) halo only contribute fully up to 
log(P 2 SAD)~ 11(8). The resolution of each simulation on 
the other hand, sets a minimum subhalo mass and thus a 
maximum value of P 2 SAD that can be trusted. These two 
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Figure 1 . Contours of the logarithm of the particle phase space 
average density ( P 2 SAD) averaged over all particles inside the 
subhaloes of four haloes of different sizes at z = 0: cluster-size 
(black), group-size (blue), Milky-Way-size (red) and dwarf-size 
(green). Two resolution levels are shown: high (solid) and low 
(dashed); see Table 1 for a summary of the simulations. The av¬ 
erage phase-space DM clustering in subhaloes is remarkably uni¬ 
versal across haloes of different masses and environments. 


scales define the innermost and outermost contours that, for 
each halo, satisfy the universality of P 2 SAD. In Fig. 1, we 
have only plotted the countours within these boundaries. 

By taking only the particles within subhaloes, conver¬ 
gence in P 2 SAD is harder to achieve due to poorer sampling 
within the volume V6, than in the case when all particles 
within the virial radius are used. This is particularly rele¬ 
vant at small scales in phase-space (log (P 2 SAD)> 11) where 
the subhaloes that dominate P 2 SAD are sampled with in¬ 
creasingly lower number of particles as P 2 SAD increases. 
Since for a fixed subhalo mass, the subhalo mass fraction 
decreases with halo mass (e.g. Gao et al. 2011), the dwarf- 
size simulation has the poorest sampling despite its better 
mass (and spatial) resolution (see Table 1). 

Fig. 2 demonstrates a different way to see the univer¬ 
sality of P 2 SAD and its natural cutoffs on large and small 
scales for a given host. It shows P 2 SAD as a function of the 
average mass enclosed by each P 2 SAD contour: 

M ave (S*) = f H(Ai',Ar')v 6 rf 3 Ax'd 3 Av', (4) 

J A(H,) 

where A(E*) is the area in (Ax, Av) defined by a fixed value 
of S*, i.e., the area under a given P 2 SAD = S* contour. 
Each halo has two effective cutoffs in this plot: one to the 
right, given by the maximum mass (size) of its subhaloes, 
and one to the left given by the minimum mass of subhaloes 
that can be resolved. In the regions that are not affected 
by these cutoffs, to a very good approximation, E oc 
which is a prediction of the spherical collapse model in a 
ACDM cosmology at small scales (Afshordi et al. 2010). 


Figure 2. The particle phase space average density, or P 2 SAD, 
in subhaloes as a function of the average mass enclosed by a 
P 2 SAD contour (see Eq. 4) in four haloes of different sizes at z = 
0: cluster-size (black), group-size (blue), Milky-Way-size (red) and 
dwarf-size (green). Two resolution levels are shown: high (solid) 
and low (dashed). 


2.1 Modelling of p 2 sad 

In Papers I and II, we used two models of P 2 SAD, one is 
simply a fitting formula where the contours of E( Ax, Av) 
are described by a family of superellipses with axes that are 
functions of E (see Eqs. 2 and 3 of Paper II): 



= 1 , 


(5) 


Previously, the functions X and V were modelled as power 
laws, but we have found that a better fit to the simulation 
data at all scales is given by the following functional form: 


X(E) = atan PoH p \V(S) = atan p 3 B P4 . 

(6) 

Notice that for log(E) ;$> P 2 (ps), these functions approach 
power laws, which are good approximations to the structure 
of P 2 SAD at small scales (see also Fig. 1 of Paper II). 

The second model is a physically motivated approach 
that combines the stable clustering hypothesis in phase 
space (Afshordi et al. 2010), the spherical collapse model, 
and tidal disruption of subhaloes (see Sec. 4 of Paper II). 
In this model, structures form according to the spherical 
collapse model with a characteristic mass, m co i, and phase 
space density, £ s (m co i), at the time of collapse: 

c _ Pchar _ 10ff(£ 3 ) . . 

a? ir ~ G 2 m col ({ s )’ U 

where the characteristic density and velocity of the collapsed 
object are (e.g. Afshordi & Cen 2002): p c h ar = 200p C rit and 
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Parametric model (Eq. 5) 

Po [Mpc/h] 

Pi 

P2 

P3 [km/s] 

P4 

P5 

P 


73.31 

-0.40 

31.04 

3.69 X 10 4 

-0.30 

7.00 

1.0 

Physical model (Eq. 10) 

-^tid 

a 

/ 

B 

K 

a 

b 3 


0.12 

1/3 

1.5 

0.192 

2.5 

0.75 

3.53 1.0 


Table 2. Best fit parameters of two different models that describe P 2 SAD\ parametric (Eq. 5) and physical (Eq. 10). We note that the 
parameter 3 is always fixed to 1.0, so these models have effectively 6 and 7 free parameters, respectively. 


cr char = er v ir = 10.HV2OO- The subhalo collapses when the 
r.m.s top-hat linear overdensity er(m co i) (mass variance) 
crosses the linear overdensity threshold 5 C ~ 1.686 at an 
epoch given by: 

( 8 ) 

The primordial phase space densities of structures are 
eventually diluted in time due to tidal stripping by a fraction 

At(ra C oi): 

S(Ak, Av) = u(m co i)€s(m co i), ( 9 ) 

which is fully determined by a tidal stripping model that has 
5 free parameters: the normalization (Atid) and slope (a) of 
the power law dependence of tidal stripping as a function of 
the ambient density, a characteristic host mass where strip¬ 
ping begins to be effective (fm co l), and the normalization 
( B ) and slope (k) that define the initial condition (pre-infall) 
as a function of the mass variance a(m co i). 

The functional form of E( Ax, Av) is given by the gen¬ 
eral solution to the collisionless Boltzmann equation under 
the stable clustering hypothesis (Afshordi et al. 2010); for 
S = const., we have: 

Ax \ ^ / Av 

a\(m co i)) \ 6 C( m coi) 

where A and £ are fully determined by the model (see Eqs. 17 
and 18 of Paper II), while a, b and /3 are free parameters that 
represent 0(1) deviations over our model to be calibrated 
to simulations. 

Eqs. 5 and 10 describe the individual contours of 
P 2 SAD from the simulation data very precisely if we fix 
/3 = 1 and fit the 6 (7) additional parameters of the paramet¬ 
ric (physical) model. Across the seven orders of magnitude of 
P 2 SAD resolved in the simulations, we can find an average 
fit that reasonably describes the structure of P 2 SAD. Ta¬ 
ble 2 shows the best fit parameters for the models (given by 
the average values of the logarithmic least squares fit to each 
contour), while Fig. 3 compares them to the simulations. 



3 CLUSTERING IN SUBHALOES AND DARK 
MATTER SIGNALS 

From the definition of P 2 SAD, we can write the real space 
two point correlation function of densities as: 

(p(x)p(x + Ax)) V{ , = J d 3 vd 3 Av (/) Ve S(Ax, Av)v 6 

= (p) Ve J d 3 Av E(Ax,Av) v 6 ,(11) 



Figure 3. Fits to the particle phase space average density 
(P 2 SAD) measured in simulations (solid lines with colors, see 
Fig. 1), using a parametric formula (magenta dashed lines, Eq. 5) 
and a physically motivated model based on stable clustering + 
spherical collapse + tidal disruption (blue dashed lines, Eq. 10). 
The best fit parameters for each model are given in Table 2. 


where (p)y 6 is the average dark matter density within the 
volume where P 2 SAD is calculated. The standard two point 
correlation function simply follows from Eq. 11: 

£(Ax)v 6 = [ d s Av E(Ax,Av) Ve - 1 , ( 12 ) 

Pm J 

where p m is the mean cosmic dark matter density. The di¬ 
mensionless power spectrum can then also be computed: 

A2( q = 2^ 3 P(fc)v 6 

To normalize ((Ai)v 6 to a cosmic volume V, we make 
the following substitution in Eq. 12: 

(pK b _ 1 My _ f su bs(V) 

Pm Pm pmV pm 

where f au bs(V) is the subhalo mass fraction within the cos¬ 
mic volume V, which is approximately given by / su b s ~ 
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k(h/Mpc) 


Figure 4. Predictions from the modelling of P 2 SAD for the 
power spectrum of DM inside subhaloes with collapse mass in 
the range 10 —6 — 2 X 10 12 Mg, normalized to the cosmic mean 
density. The lower value represents the canonical minimum self¬ 
bound mass for 100 GeV WIMPs, while the maximum mass 
is given by the minimum value of P 2 SAD that we can trust 
(log(P 2 5j4Z)) m i n = 6, see Fig. 3). We also show the halo model 
predictions for the power spectrum in haloes in the same mass 
range using the mass function from Schneider et al. (2013), and 
either a NFW profile with two different concentration-mass rela¬ 
tions (SCP14, Sanchez-Conde & Prada 2014, dashed black, and 
OA15, Okoli & Afshordi 2015, dashed blue), or a Einasto profile 
with parameters given in Klypin et al. (2014) (K14, dashed red). 

/subs(14iw)/h ~ 0.12, where fh is the halo mass fraction, and 
fsuhs (Vmw ) is the subhalo mass fraction in a Milky-Way-size 
halo. For m m i n = 10 _6 Mg, fh ~ 0.8 (e.g. Angulo & White 
2010) and / su bs ~ 0.15 (e.g. Springel et al. 2008). 

In Fig. 4 we show the dimensionless power spectrum in 
subhaloes as predicted by P 2 SAD. To calculate the power 
spectrum, we have used the FFTL0G code (Hamilton 2000) 
to compute the Fourier transform of ((Ar)y 6 . To make 
this plot we have used our models of P 2 SAD and compute 
((Ai)v 6 using as limits for the integral in velocities those 
that correspond to the values of S. The minimum is given 
by the level at which we can trust the models in its com¬ 
parison with simulations (log^^ADJnm, = 6, see Fig. 3), 
while the maximum is given by the CDM particle model one 
wishes to study. As an illustration, we have chosen a 100 
GeV WIMP with a canonical value for the minimum self- 
gravitating mass for collapse m m i n = 10 6 M 0 . Since our 
physical model for P 2 SAD allows us to connect a collapse 
mass with a value of H (see Eq. 9) , we can easily interchange 
m co i with Av (for a given Ax) in Eq. 12. In this way, Fig. 4 
shows the DM clustering in subhaloes with collapse masses 
in the range 10 -6 — 2 x 10 12 Mq. Remarkably, both models 
of P 2 SAD give a very similar result, despite the fact that 
they are only calibrated within the scales resolved in the 
simulation. At large scales the clustering is artificially sup¬ 


pressed by the minimum limit we imposed on P 2 SAD, while 
at small scales, A 2 (fc) ~ const., due to the physical cut-off 
given by m min . 

In Fig. 4, we also show the DM power spectrum ac¬ 
cording to the halo model (e.g. Seljak 2000), for halo masses 
in the same range as the one used for P 2 SAD : 1CK 6 — 2 x 
10 12 Mq. For the ingredients of the halo model we used the 
halo mass function for a Planck cosmology using the code by 
Schneider et al. (2013), and two different halo profiles, NFW 
profile with a shallow concentration mass relation (black 
dashed line, Sanchez-Conde & Prada 2014) and Einasto pro¬ 
file with an effective steep concentration mass relation (red 
dashed line, Klypin et al. 2014); we also show the result for 
NFW profile with the concentration mass relation given by 
Okoli & Afshordi (2015) (blue dashed line). It is important 
to remember that since P 2 SAD was defined only with par¬ 
ticles inside subhaloes, our results cannot be directly com¬ 
pared with the halo model, which is defined by the DM clus¬ 
tering inside (and across) haloes. At small scales however, 
the clustering is dominated by the so-called one-halo term, 
which is defined as correlations between particle pairs within 
collapsed haloes (subhaloes). Because of this, P 2 SAD pre¬ 
dicts the dominant contribution to the small-scale A 2 (fc), 
or £(Aa;), for all matter, albeit with a normalization that 
depends on the uncertain total DM mass contained in sub¬ 
haloes (see Eq. 14). 

Our prediction of the non-linear power spectrum at 
small scales significantly decreases the uncertainty us¬ 
ing other methods/extrapolations (see e.g., Fig. 3 of 
Sefusatti et al. 2014 for a plot containing several of these 
predictions). In the halo model for instance, this uncertainty 
is ultimately connected with the sensitivity of the extrap¬ 
olation (towards low-masses) of the concentration-mass re¬ 
lation on the assumptions regarding the structure of haloes 
in the resolved regime. This is illustrated in the difference 
between the cases SCP14 (black dashed line) and K14 (red 
dashed line) in Fig. 4. These two cases are based on calibra¬ 
tions with similar simulation data but using different pro¬ 
files (NFW and Einasto, respectively). Their predictions for 
the concentration-mass relation at low masses, and therefore 
on the power spectrum, are quite different. Our method on 
the other hand, is based on the modelling of a 2D function 
( P 2 SAD ), which is resolved over many orders of magnitude. 
The fact that our two different models (physical and para¬ 
metric) give a very similar result for the power spectrum at 
small scales (magenta and blue solid lines in Fig. 4) indi¬ 
cates, in contrast to the halo model , a low sensitivity to the 
way the fit (calibration) is done in the resolved regime. We 
argue that this is an advantage of our method compared to 
the halo model. 

In this work we have concentrated in analysing DM-only 
simulations. In the presence of baryons, the DM phase space 
clustering is expected to be modified as the DM responds dy¬ 
namically to the assembly of galaxies. We anticipate that the 
structure of P 2 SAD will be modified mainly at the scales 
where subhaloes contain condensed baryons (gas and stars). 
The focus of this work is however at the scales of low-mass 
subhaloes, which are devoid of cold condensed gas and stars. 
Therefore, baryons can only impact P 2 SAD at these small 
scales indirectly through the orbital interaction of subhaloes 
with the host galaxy and its satellites. This is mostly im¬ 
portant for subhaloes that at any point in their orbit were 
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Figure 5. Predicted global subhalo boost (i.e. within r 2 oo) to the 
annihilation rate of a host halo as a function of its mass, according 
to our physically motivated model of P 2 SAD (Eq. 10). The two 
cases are for two different minimum subhalo masses as shown in 
the legend. The annihilation rate in the smooth host halo is com¬ 
puted assuming it follows a NFW profile with the concentration- 
mass relation given in Sanchez-Conde & Prada (2014), with a 
scatter of 0.14 dex. The star symbol is the boost for the Aquarius 
A halo, which has a concentration of ~ 19. 


near the centre of the host galaxy. Since in this work we 
average P 2 SAD globally within the virial radius of the host 
halo, the average is dominated by the abundant subhaloes 
near the virial radius, most of which have recent infall times. 
Thus, at small scales, P 2 SAD (averaged in this way) is not 
likely to be affected by the presence of baryons. 


3.1 Global subhalo boost to DM annihilation 

The DM annihilation rate R ann (number of self-annihilation 
events per unit time) in a region of physical volume V (phase 
space volume Ve) is given by the limit of P 2 SAD on small 
distances. The global subhalo boost to the annihilation rate 
of the smooth host halo with a density profile p(x) is (Eqs. 
4, 5 and 9 of Paper II): 

Kl h n _ Mv 6 f d 3 Av(cr ann v) limArE— >0 5(As, At>)y 6 
-R|nn° th fy d 3 xp 2 (x ) (a ann v) 


d subs _ onn n A/f (^ ann *•') 

R&nn — g£3 -00p C rit,oMv 6 ^ m 2 

/* m max 

x / u( m co\)m~£d(m 2 ol a 3 (m co i)),{ 16 ) 

J m min 

where m x is the mass of the dark matter particles, p C rit is the 
critical density, and m min ( max ) is the minimum (maximum) 
collapse mass for the subhaloes contributing to P 2 SAD. 

Using the parameters in Table 2, we can use Eqs. 15 
and 16 to compute the global subhalo boost to the annihi¬ 
lation rate of host haloes of any mass (Fig. 5). To compute 
the total mass within subhaloes, M\> e , we use the formula 
given in Gao et al. (2011). We show two values of ra m i n , and 
assume a NFW profile for the host halo with a concentration- 
mass relation as given Sanchez-Conde & Prada (2014). No¬ 
tice that the boost strongly depends on the concentration of 
the smooth host halo, e.g., if we take the Aq-A simulation, 
which has a very high concentration for its mass, we get a 
lower value for the boost (red star symbol), consistent with 
our previous results (see Fig. 5 of Paper II). 

Our predictions for the subhalo boost are fairly 
high compared to, for example, those made by 
Sanchez-Conde & Prada (2014) (see their Fig. 2). The 
reason for this is apparent in Fig. 4, where we show the 
wide range of predictions for the power spectrum at small 
scales (and thus indirectly the subhalo boost) according to 
the halo model. This uncertainty is ultimately connected 
with the uncertainty on the concentration-mass relation 
at small masses. The work by Sanchez-Conde & Prada 
(2014) is in the low end of these predictions and thus, 
predicts a lower subhalo boost than our model. It is 
important to mention that the latter work assumes that the 
concentration of a subhalo is the same as that of an isolated 
halo of the same mass, while simulations have shown that 
subhaloes are typically more concentrated. For instance, 
in a Milky-Way-size halo, subhaloes are on average ~ 30% 
more concentrated (see e.g. Fig. 28 of Springel et al. 2008). 
Given the strong dependence of the annihilation rate on 
the subhalo concentration, the subhalo boost is expected 
to be a factor of 2 — 3 higher than the estimate given by 
Sanchez-Conde & Prada (2014) once this effect is taken 
into account (see e.g. Bartels & Ando 2015, for a recent 
detailed calculation of this effect). This would bring their 
prediction much closer to ours. 


4 CONCLUSIONS 

The particle phase space average density ( P 2 SAD ) is a 
coarse-grained phase space density that contains a wealth 
of information on the clustering of matter (Papers I and II). 
By studying DM structures with this novel statistics, we 
have found a new universality in the gravitational clustering 
of DM within subhaloes: 


where (<7 ann u) is the product of the annihilation cross section 
and the relative velocity between pairs, and ( a ann v ) is its 
average over the velocity distribution of DM particles. Using 
the physical modelling of P 2 SAD, we can transform the 
numerator in Eq. 15 into an integral over m co i, which in the 
case (cr a nnu) = const., takes a simple form (Eq. 32 of Paper 


• The structure of P 2 SAD averaged over particles inside 
subhaloes is universal across haloes in a wide range of masses 
covering 3 orders of magnitude in physical and velocity sep¬ 
arations, from the subhaloes of dwarf-size haloes to those of 
cluster-size haloes (see Fig. 1). 

• The functional form of P 2 SAD (averaged in subhaloes) 
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can be described by a family of superellipse contours in the 
(Aa;, Av) plane with axes that can be modelled by simple 
parametric formulae (see Eq. 6), or by a physically motivated 
model based on the stable clustering hypothesis in phase 
space, the spherical collapse model, and tidal disruption of 
subhaloes (see Eq. 10). 

We have exploited this new universality of DM clustering to 
accurately predict the highly non-linear DM power spectrum 
all the way down to the decoupling limit of CDM particles 
(Fig. 4). This prediction can be used to study several poten¬ 
tial DM signals, e.g., the subhalo boost to the annihilation 
of DM for different host halo masses. A code to compute 
P 2 SAD using our physical model is available upon request. 
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